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Abstract 

Frequency dependent selection and demographic fluctuations 
play important roles in evolutionary and ecological processes. 
Under frequency dependent selection, the average fitness of the 
population may increase or decrease based on interactions be¬ 
tween individuals within the population. This should be reflected 
in fluctuations of the population size even in constant environ¬ 
ments. Here, we propose a stochastic model, which naturally 
combines these two evolutionary ingredients by assuming fre¬ 
quency dependent competition between different types in an 
individual-based model. In contrast to previous game theoretic 
models, the carrying capacity of the population and thus the pop¬ 
ulation size is determined by pairwise competition of individuals 
mediated by evolutionary games and demographic stochastic- 
ity. In the limit of infinite population size, the averaged stochas¬ 
tic dynamics is captured by the deterministic competitive Lotka- 
Volterra equations. In small populations, demographic stochas- 
ticity may instead lead to the extinction of the entire population. 
As the population size is driven by the fitness in evolutionary 
games, a population of cooperators is less prone to go extinct 
than a population of defectors, whereas in the usual systems of 
fixed size, the population would thrive regardless of its average 
payoff. 


Introduction 

All natural populations are composed of a finite number of in¬ 
dividuals. These individuals can reproduce, interact, die or mi¬ 
grate, which leads to changes in the population size over time. 
In many theoretical models, it is convenient and possible to ne¬ 
glect the effect of demographic fluctuations by assuming infinite 
populations when population sizes are sufficiently large (1) or by 
assuming constant population size as in the Moran or Wright- 
Fisher process (2). However, such simplifications may be invalid 
when considering additional ecological processes including os¬ 
cillations in population size of predator and prey systems (3- 
5), periodic fluctuations and outbreaks of infectious diseases in 
humans (6), or chaotic dynamics under multi-species interac¬ 
tions (7). The Lotka-Volterra equations provide a deterministic 
description of the abundances of species as continuous densi¬ 
ties but they are not designed to include the impact of random 
drift. Theoretical models coupling changing population size and 
stochastic dynamics arising from individual based models have 
become more popular only recently (8-12). 

The Lotka-Volterra equations naturally take frequency depen¬ 
dent selection into consideration. Under frequency dependent 
selection, the fitness of a given type (or species) depends on 


the composition of the entire population (or community) (13). 
Different kinds of frequency dependence can lead to different 
dynamical patterns. Most prominently, negative frequency de¬ 
pendent selection can result in a stable coexistence of differ¬ 
ent types (14-18). One elegant way to describe such frequency 
dependence is through evolutionary game dynamics (13, 19). 
Evolutionary game theory has extensively developed the theory 
of stochastic dynamics in the past decade (2). However, most 
progress has been accomplished for constant, finite population 
sizes, which is mathematically convenient, but not always re¬ 
flects biologically appropriate scenarios. Game theoretic models 
that take changing population sizes into account mostly focus 
on deterministic dynamics (20) - similar to traditional ecologi¬ 
cal models. Here we introduce a simple and elegant model of 
stochastic evolutionary game dynamics that explicitly allows for 
changing population size through a natural interpretation of pay¬ 
offs in terms of competition between individuals. 

Traditional game theoretic models assume that individuals ob¬ 
tain payoffs from interactions with other members in the popu¬ 
lation, which are translated into fitness. Individuals with higher 
fitness are assumed to have more offspring and hence repro¬ 
duction is frequency dependent. Conversely, death rates are as¬ 
sumed to be constant or normalized such as too keep the popu¬ 
lation size constant. In contrast, here we focus on a microscopic 
description of the dynamics in terms of reaction kinetics equa¬ 
tions. Assuming constant birth rates and frequency dependent 
death rates allows to interpret payoffs directly and naturally in 
terms of reaction rates, where selection acts on survival rather 
than reproduction. The present setup lends itself to a straight 
forward derivation of the deterministic dynamics in the limit of 
large population sizes in the form of competitive Lotka-Volterra 
equations, but equally allows to model the stochastic dynamics 
in finite populations of variable sizes, which may even lead to the 
extinction of the population. 


Model and Results 
Stochastic dynamics 

Most models for stochastic evolutionary game dynamics con¬ 
sider a fixed population size, such that every birth is balanced 
by the death of another individual (2). Simply decoupling birth 
and death events in such models leads to random fluctuations 
in the population size and thus eventually to stochastic extinc¬ 
tion (21). Instead, here we propose a framework based on the 
microscopic processes of birth, death and competition. For sim¬ 
plicity, we focus on two types of individuals, X and Y, but note 
that the generalization to arbitrary numbers is straight forward. 
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Every individual reproduces 


+ X and Y^V + Y (1) 

at constant rates and Ay^yy, respectively, and dies 

X ^ 0 and y ^ 0 (2) 

also at constant rates and A^^o, respectively. Competitive 
interactions result in four more processes 

x + x^x x + y^y [3a] 

x + y^x y + y^y. [3b] 

In the simplest case, all competition rates are equal, such that 
two randomly chosen individuals compete for survival (9). 


Competitive seiection 


The most natural way to introduce evolutionary games in the 
above framework is to relate the competition rates in Eq. (3a) to 
a payoff matrix 
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which determines the strength of competition between two X, Y 
individuals as a, b, c, and d such that individuals obtaining higher 
payoffs are less likely to die in competitive interactions. More 
specifically, we assume that reaction rates are the inverse pay¬ 
offs scaled by M, a unit for controlling population size. 


\ — ^ \ — 

y\xx^x — '^xy^y — , 

Axy^x — '^yy^y — dM' 

Thus, if, for example, an X and a Y individual compete, the X 
dies with a probability proportional to l/(6M) and the Y pro¬ 
portional to l/(cM). This requires that a,b,c,d > 0 in order 
to remain meaningful in terms of reaction rates. In traditional 
models, payoffs are associated with reproduction, whereas here 
they refer to the probability of surviving competitive interactions. 
In both scenarios, high payoffs result in increased reproduc¬ 
tive output over the lifetime of an individual. Naturally, reaction 
rates could follow a different functional dependence, for example 
Axx^x = exp(-aM), which would lead to qualitatively similar re¬ 
sults as long as rates decrease with increasing payoffs, but with¬ 
out the restriction of positive payoffs. The scaling term M deter¬ 
mines the frequency of competitive interactions as compared to 
birth or death events, Eqs. (1) & (2). As long as population sizes 
are much smaller than M, competition is rare and most events 
are births or deaths. In large populations, however, competition 
becomes common and results in density dependent regulation 
of the population size. 

Implementing evolutionary games through competition is, of 
course, just one approach to link payoff matrices to reaction 
rates. Intuitively, it is tempting to assume that evolutionary 
games determine the birth rates instead because payoffs then 
more directly reflect fitness advantages. However, this requires 
microscopic reactions of the form X-i-X^-X-i-X-i-X, which 
occurs at rate a, etc. Such interactions seem more appealing 
in sexually reproducing populations rather than for the more tra¬ 
ditional models, which focus on one sex only or are based on 
asexual reproduction. More importantly, however, even when in¬ 
cluding competition at fixed rates, this setup remains inherently 
biologically not meaningful, because it either results in extinction 
or indefinite growth of homogeneous populations (see Support¬ 
ing Information, SI, Sec. 1). In natural systems, there can be 
positive feedbacks between different types (22), but they typi¬ 


cally refer to different systems where other effects, such as pre¬ 
dation, guarantee that the population size remains finite. Here, 
we focus on a competitive system with negative feedback in¬ 
stead. 


Large population size 


The reaction based system above can be formulated in terms 
of a continuous-time Master equation, see SI, Sec. 2. For large 
M, a measure of the equiibrium densities, this equation can be 
approximated by a Fokker-Planck equation, which describes the 
dynamics of the probability distribution of the system (23). When 
the population densities approach the equilibria, we recover de¬ 
terministic rate equations from the microscopic processes de¬ 
fined in Eqs. (1)-(3a) 
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where x and y denote the density of individuals of type X, Y and 
ra, = Ax^xx - Ax^o, Xy = Ay^yy “ Ay^Q indicate the intrinsic 
growth, i.e. the net growth rates from birth and death events, 
Eqs. (1) & (2). Note that the deterministic limit can be derived 
directly based on the law of mass action. Even though only the 
net growth rates, Xx and Xy, enter Eq. (5), it is important that the 
stochastic description does not lump the two processes together 
in either decreased net birth or net death rates. In particular, if 
spontaneous death events, Eq. (2), are dropped (or absorbed in 
reduced birth rates, Eq. (1)) such that deaths occur only due to 
competition, Eq. (3a), then populations would never go extinct in 
the stochastic formulation because the last individual standing 
would remain immortal. 

In order to recover the familiar form of the competitive Lotka- 
Volterra dynamics (24, 25), we factor out Xx in Eq. (5) and set 
0=1 (without loss of generality), 
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Here K = XxM simply denotes the carrying capacity of X types 
and XydM the corresponding carrying capacity of Y types. In 
the absence of Y types Eq. (6a) reduces to the logistic equation 


X 



[7] 


which forms the basis for r-X-selection theory (26), where the 
carrying capacity K is independent of the intrinsic growth rate, 
Xx- However, according to Eq. 5, K is an emergent quantity 
determined by the population’s environment (27-30), which cru¬ 
cially includes all members of the population together with their 
ecological interactions. For example, if a mutant type Y dou¬ 
bles its intrinsic rate of reproduction as compared to the resident 
X, Xy = 2xx, then the mutant type readily displaces the resi¬ 
dent and reaches its carrying capacity at twice the density of the 
resident, = Xy/xxK = 2K, assuming that all other environ¬ 
mental parameters remain the same. This conclusion does not 
only follow from the microscopic description of relevant biological 
processes but has also been observed in experimental settings 
(27, 28, 31, 32). Moreover, in the following we show that the 
notion of a carrying capacity becomes even more challenging in 
populations of multiple types. Thus, we use the density of indi¬ 
viduals at equilibrium instead of carrying capacity in heteroge¬ 
nous populations, where the total equilibrium density is Kcox and 
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the densities of type X and type Y individuals at this equilibrium 
are and respectively . 

If < 0 or Ty < 0, the corresponding type will invariably de¬ 
cline and disappear. Including competition only speeds up their 
demise. For example, this applies to the predators in the fa¬ 
mous, oscillating Lotka-Volterra predator-prey dynamics. Since 
negative interaction rates are not meaningful, at least three dif¬ 
ferent types are required to observe oscillations in a competitive 
system (25, 33). 

Equilibria of the deterministic system 

The deterministic mean-field dynamics of our model serves as a 
valuable reference for the underlying stochastic evolutionary pro¬ 
cess. Birth and death rates, Eqs. (1) & (2), may differ for different 
types but for the sake of simplicity and to highlight effects arising 
from evolutionary games, we discuss interactions of two types 
of individuals, X and Y, with = Vy = r. The evolutionary fate 
of each type depends on a combination of the strength of intra¬ 
type competition (a and d) and inter-type competition {b and c). 
In general, we can classify three different scenarios based on 
the payoff matrix: (i) For a > c and b > d, type X individuals in¬ 
variably achieve higher payoffs (i.e. longer life expectancy) than 
Y types and hence type X dominates type Y. Similarly, type Y 
dominates type X whenever a < c and b < d, see Fig. la. The 
prisoner’s dilemma is the most prominent example of a domi¬ 
nance game (19). (ii) For a > c and b < d both types are at 
a disadvantage compared to the other type when rare. This re¬ 
flects coordination games such as the stag-hunt game (34). (iii) 
Finally, for a < c and b> d both types have an advantage when 
rare but are at a disadvantage when abundant. Thus, an interior 
equilibrium exists where the two types co-exist, see Fig. 1b. The 
hawk-dove or snowdrift games are examples of such scenarios 
(35). Note that even though the classification of the dynamics 
for two types is based on their payoffs in the same way as in 
the classical replicator dynamics (and the stability remains the 
same, see SI, Sec. 3), the position of the rest points in our de¬ 
terministic system are naturally different. For example, in the 
replicator dynamics a coexistence game as in Fig. 1b exhibits 
a stable rest point at x* = {d- b)/{a - b - c + d) = 1/2. In 
contrast, according to Eq. (6), the frequency of X at equilibrium 
is X* = 10/13. An intuitive reason for this increase in the rel¬ 
ative abundance of X is, that the total number of individuals at 
the mixed equilibrium is lower than the carrying capacity for a 
population of only X types. 

Therefore, the only possible equilibria are either homogenous 
X or X populations or a stable heterogenous mixture of the two. 
According to Eq. (5) with r^, = = r, the densities of indi¬ 

viduals at the three equilibria are AT"' = aMr and = dMr 
as well as Xcox = K^ox + K^ox with K^^x = ^nd 

Xlox = ^1’’ which can be rewritten as Kcox = K"" + 

Note that in the co-existence equilibrium the den- 
sity of individuals of each type is always lower than its carrying 
capacity in isolation, i.e. K^ox < and < K^. Flowever, 
the total number of individuals in mixed equilibria can either ex¬ 
ceed or fall short of homogenous carrying capacities: if 6 > a, 
Kcox > AT"; if & < a, ATcox < if c> d, Koox > K^', if c < d, 
Kcox < K^- More specifically, in co-existence games b > a and 
c > d holds such that the total number of individuals is highest in 
the mixed equilibrium, Kcox > K'^,K^. Conversely, the reverse 
ranking is impossible: it would require b < a and c < d but this 
refers to coordination games where the mixed state is unstable 
and the population approaches either one of the homogenous 
configurations. Of course, in the corresponding stochastic re¬ 
alisations the population size at equilibrium is not fixed and in¬ 
stead fluctuates around the carrying capacity, see Fig. 1. For 
identical birth and death rates, the evolutionary game controls 


the relative growth or decline of the two types through competi¬ 
tion, but regardless of the game the numbers of both types can 
increase or decrease if the current state of the population is far 
from equilibrium. All possible rankings of equilibrium densities 
are summarized in Table 1 in the SI. 

Stochastic simulations 

In contrast to the deterministic equilibrium predictions, the only 
evolutionary outcome in stochastic simulations is the eventual 
extinction of the entire population - all other states are tran¬ 
sient. Fortunately, the expected times to extinction rapidly grow 
with the density of individuals in equilibrium, controlled by M, 
Fig. 2. Flence, predictions based on deterministic dynamics, see 
Eq. (5), keep providing valuable insights for the stochastic dy¬ 
namics, see Eqs. (1) - (3a), especially for large population sizes. 
Substantial quantitative and even qualitative differences arise as 
illustrated Fig. 1 for three characteristic types of interactions. The 
stochastic dynamics is implemented through the Gillespie algo¬ 
rithm (36). For the relatively large carrying capacity X’' = 1000 
in Fig. 1, each realisation of the stochastic dynamics fluctuates 
around the deterministic trajectory. Fluctuations represent an 
integral part of natural populations and hence stochastic evolu¬ 
tionary models provide a more natural way to study evolutionary 
trajectories, especially to capture the interplay between ecologi¬ 
cal and evolutionary processes (10, 11,37, 38). 

In small populations, competition for survival is weak and the 
dynamics is mainly determined by the intrinsic growth rate, r, i.e. 
individual birth and death events. Consequently, small r results 
in higher stochasticity but also tends to decrease the number of 
individuals at equilibra, which further amplifies the effect. As the 
population grows and approaches its carrying capacity, competi¬ 
tion becomes increasingly important and competition rates, see 
payoff matrix Eq. (4), also control the size of fluctuations. Strong 
competition (small payoffs) reduces stochasticity, but also tends 
to decrease the population size, which may offset the reduction 
in terms of fluctuations. 

Furthermore, if the numbers of the two types in the stochas¬ 
tic process are far away from the deterministic equilibrium, even 
the averaged stochastic dynamics can be very different from de¬ 
terministic predictions. For example, in the deterministic case 
a dominant mutant always succeeds in invading and eventually 
fixating in the population. In contrast, in the stochastic case a 
single mutant often fails to invade and fixate even if it is domi¬ 
nant. Note that fixation of a single mutant becomes even less 
likely in larger populations - in spite of the fact that fluctuations 
decrease and the deterministic dynamics is recovered in the limit 
of large populations. 

Extinction 

In eoological models the risk of extinction of a population due 
to demographic stochasticity has recently received considerable 
attention, see e.g. (10, 37, 38). Evolutionary game theory mod¬ 
els demonstrated that stochastio fluotuations are important de¬ 
terminants for the fixation and extinction of individual traits even 
under constant population sizes (2) but remained unable to ad¬ 
dress the more dramatic possibility of the extinction of the entire 
population. 

In a dominance game, the deterministic dynamics predicts 
that the dominant type invariably takes over the entire popula¬ 
tion and approaohes its carrying capacity, see Fig. 2a. Flowever, 
in stochastic models the two types X and Y may go extinot 
sequentially due to fluctuations. Especially when starting from 
small populations, the extinction probability is not negligible, see 
Fig. 2b. Note that the (cumulative) extinction probability oon- 
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Figure 1: The stochastic dynamics fluctuates around the deterministic predictions in large populations (black lines). Each panel 
refers to one of the three characteristic classes of interactions as determined by the ranking of payoffs. Evolutionary 
trajectories for X (blue) and Y (red) types are shown for initially small (xo = 10, yo = 90, pale colours) and large 
(aio = 1000, yo = 1000 saturated colours) populations. Regardless of the game both types tend to increase for low initial 
population sizes, whereas both types decrease in densities at large initial sizes. In all cases the deterministic predictions 
agrees well with the stochastic dynamics in that the stochastic trajectories fluctuate around the deterministic average 
(parameters M = 2000, = 0.6, K-to = Ay^o = 0.1, which translates into = 1000 and = 250 in 

(a) and = 500 in (b), (c)). 



Time 


Figure 2: Stochastic dynamics of a dominance game in small 
populations, (a) Three realisations of the stochastic 
dynamics in a small population with two types, X (blue) 
and Y (red) for M = 1000. Stochastic trajectories 
deviate significantly from the deterministic dynamics 
(solid black lines). For example, in the third realisa¬ 
tion both types go extinct, first the dominant type AT 
followed type Y a little later, (b) Cumulative probabil¬ 
ity for the extinction of the entire population over time, 
averaged over 10® realisations for different M (param¬ 
eters Aa:— yxx — Ay — yyy — 0.6, \x —^0 ~ Ay —^0 ~ 0.5, 
xo = 1,2/0 = 9. This yields TsT* = 100 and = 25 in 
the upper panel). 


verges to one as time goes on regardless of the population size. 
Only for larger populations extinction typically takes much longer 
and the extinction probability increases slowly with time. The 
smaller the population size - due to small carrying capacities, 
fluctuations, or initial configurations - the higher the risk of ex¬ 
tinction. 

In Fig. 2a, the dominant X type has a higher carrying capacity 
than the Y type and hence, provided that the X mutant success¬ 
fully invades and fixates, the population ends up more persistent. 
In the third realisation in Fig. 2a, the X type is lost first and the 
entire population vanishes soon after. An essential factor for the 
persistence of small populations are the birth and death rates, 
Eqs. (1) & (2), or, more specifically, their difference: for large in¬ 
trinsic growth rates, small populations are more likely to escape 
extinction because the population can more readily recover and 
return to its carrying capacity. 

In co-existence games, the two types X and Y typically co¬ 
exist and their densities fluctuate around their respective (deter¬ 
ministic) densities of individuals in equilibrium, X5,x and 
see Fig. 1b. Whenever stochastic fluctuations drive one type to 
extinction, the carrying capacity of the remaining type changes 
to or K^, respectively. Thus, the density of individuals of 
the entire population changes and is driven by the interplay of 
demographic fluctuations arising from intrinsic growth rates and 
the strength of competition, see Fig. 3. 

In this case the extinction of the entire population is a two- 
step process: first one types goes extinct - typically the type 
with the lower density of individual in equilibrium - and then the 
population fluctuates around the homogenous carrying capacity 
of the remaining type such that the extinction dynamics is now 
governed by the stochastic equivalent for the logistic growth of a 
single type, which is well understood (10, 24, 37). 

Snowdrift game 

As a concrete example of a co-existence game, we consider the 
stochastic dynamics of the snowdrift game (35) (or, equivalently, 
the hawk-dove game, (19)). In the snowdrift game two individ¬ 
uals need to finish a task, which provides benefits I3 to both. 
The costs of the task, 7, are shared equally if both cooperate, 
i.e. participate in completing the task. If only one participates, 
the cooperator has to bear the entire costs but the defector still 
receives the benefits. Finally, if both defect and refuse to partic- 
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Figure 3: Stochastic dynamics and large demographic fluctu¬ 
ations in a co-existence game. Top (M = 2000): 
According to the deterministic dynamics (solid black 
lines) the two types X (blue) and Y (red) co-exist at 
their respective equilibrium, K^ox and K^ox- The three 
stochastic realisations, however, reveal a different pic¬ 
ture: In the first case, type X goes extinct first, which 
benefits type Y and results in a higher carrying ca¬ 
pacity than the density of type Y individuals in the co¬ 
existence equilibrium, K’^ > K^^x- In the second case, 
type Y goes extinct first and the new carrying capac¬ 
ity iT* even exceeds the whole population size in the 
co-existence equilibrium > Kcox- In the third case 
(inset), both types go extinct early on. Pale lines repre¬ 
sent the simulated standard deviation (averaged over 
500 independent simulations with > 10® data points 
each). Bottom: The standard deviation in the coex¬ 
istence equilibrium increases with M, but differently 
for the two types (parameters = >^y^yy = 0.6, 

Aa,^o = Ay^o = 0.5, Xf) = 10, j/o = 10). 


ipate, the payoffs for both individuals are close to zero. Hence, 
the payoffs of cooperators, AT, and defectors, Y, are a = 
b = ,0-7, c = / 3 , and d = e, where ,0-7 > e > 0 and ^ >'r > 0. 
Note that in the limit e 0 the death rate of Y types due to com¬ 
petition diverges and they are no longer able to persist in isola¬ 
tion. In order to maximize its own gain, each player should do 
the opposite of what its opponent does. A population of cooper¬ 
ators fluctuates around M. A defecting mutant 

has a selective advantage and hence is likely to successfully in¬ 
vade and the population typically starts fluctuating around a new 
equilibrium, where the total density of individuals is lower- from 
b < a follows Kcox < K"^, see Table 1 in the SI. For sufficiently 
small £, the conditions for the ranking A'”’ > K^ox >Ky > K^cox 
are satisfied. Thus, even though defecting mutants are favoured, 
their abundances in the co-existence equilibrium is even lower 
than when in isolation. Consequently, stochastic fluctuations 
are more likely to eliminate defectors and re-establish cooper¬ 
ation at the original carrying capacity, For example, for 
M = 2000, = Ty = 0.5, ;0 = 1.5,7 = l,e = 0.05, the den¬ 

sity of cooperators K^ox ~ 964 in the mixed equilibrium is much 
larger than that of defectors ATcox ~ 18 and hence the odds of 
persistence are clearly in favour of cooperators. 

Discussion 

Demographic fluctuations based on ecological interactions cap¬ 
ture important aspects and represent crucial determinants 
of evolutionary trajectories, especially in smaller populations. 
Here, we introduce a simple stochastic framework built on the 
microscopic events of birth, death and competition. This frame¬ 
work admits a simple yet elegant way to implement evolutionary 
games through payoff based competition rates, which results in 
selection on survival instead of the more traditional fecundity 
based selection (30). This yields a stochastic model for evo¬ 
lutionary games in populations of changing and fluctuating finite 
size. In the limit of infinite population sizes, this framework recov¬ 
ers the deterministic dynamics of the competitive Lotka-Volterra 
equations and hence allows to pinpoint and emphasize differ¬ 
ences that arise due to stochastic effects. The deterministic limit 
of the stochastic framework also highlights that, in contrast to the 
classical r-AT-selection theory, the ecological carrying capacity 
of a population is an emergent quantity (27, 28), which depends 
on the population configuration and is determined by the under¬ 
lying processes of birth, death and competition. In particular, 
mutations that alter the rates of any of these processes trigger 
a change in the (deterministic) carrying capacities of the mu¬ 
tant population, provided that it succeeds to take over, or of the 
mixed population in the case of co-existence. Our model implies 
that adaptation is not a simple process of accumulating benefi¬ 
cial mutations with higher carrying capacities in isolation, but in¬ 
stead an adaptive process that can favour invasion and fixation 
of mutations that are disadvantageous for the entire population, 
e.g. evolutionary suicide(39). 

Dominant mutations are bound to take over with certainty un¬ 
der deterministic dynamics. However, in the stochastic scenario, 
the chances for a single beneficial mutant to successfully invade 
and take over remains small, even for a dominant strategy. At 
first it might be surprising that the chances of success decrease 
for increasing population sizes - despite the fact that the limit 
of large populations recovers the deterministic dynamics. But, 
of course, in this limit the mutant density converges to zero, 
which resolves the apparent contradiction. Similar results can 
be found in classical models of finite populations with constant 
fitness values (40). Here we investigated stochastic dynamics 
in well-mixed populations but a natural extension is to consider 
spatial dimensions, which may increase stochastic effects due 
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to small local sub-populations (41-43). 

Here, we have focussed on the paradigmatic case of one pop¬ 
ulation and two types, but It Is straight forward to extend the 
framework to Include multiple types. In a population with three 
types, oscillations can persist in the stochastic process, while the 
deterministic limit suggests cycles spiraling towards an internal 
equilibrium (44, 45), see Fig.1 in the SI. Moreover, our frame¬ 
work easily extends to group interactions, such as public goods 
games, by allowing for competitive interactions that involve more 
than two individuals. However, in either case, the number of mi¬ 
croscopic interactions tends to increase rapidly and hence ham¬ 
pers a more general yet compact and intuitive presentation. 

The stochastic framework also emphasizes that in the long 
run populations invariably go extinct, which means that the de¬ 
terministic equilibria merely indicate fleeting states - albeit the 
expected time to extinction can be exceedingly long, especially 
for larger populations. Therefore, it remains reasonable to con¬ 
sider the deterministic predictions as a baseline superimposed 
by fluctuations of stochastic realisations. At the same time it is 
crucial, especially in smaller populations, to consider the persis¬ 
tence of individual traits or the viability of the entire population. 
For example, in the snowdrift game an invasion attempt by de¬ 
fectors triggers ecological feedback, which alters the carrying 
capacities in favour of cooperators, such that stochastic fluctu¬ 
ations help to eliminate defectors and re-establish homogenous 
cooperation. 
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1 Relating payoffs to birth 

If payoffs affect birth instead of competition, the corresponding reactions can be written as 

x+x Ax+x+x 

X+Y \X+Y+X 
x+Y Ax+y+y 
y + y 4y + y + x. 

Using the same principle as in the main text, i.e. individuals with higher payoffs have advantages 
in pairwise interactions, the birth rates should be an increasing function of the payoff elements 
a, b, c, and d. The simplest choice is that birth rates equal payoffs. We assume constant intrinsic 
death rate Xd for both types, which implies 

X ^ 0 and y ^ 0. 

Following a logistic growth model, we assume neutral competition 

X + X ^ X 
x + y ^ X 
x + y ^ y 
y + y ^ y. 
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Combining the ten reactions above, we obtain the deterministic rate equations 

X = x(ax + by- Ac(x + y) - Ad), 

y = y{cx + dy-\c{x + y)-\d), ( 1 ) 

where x and y are the number of individuals of type X and type Y. We denote the frequencies 
by n = x/{x + y) and v = y/{x + y). With Eqs. (1), we obtain the change of the frequency of 
type X as 

xy — yx 

u = --^ 

(x + yY 

= {x + y) u{l — u)[{a — h — c + d)u + b — d]. (2) 

N Replicator dynamics 

The dynamics is equivalent to the replicator dynamics scaled by the total population size N. 
As N changes over time, this can be considered as a dynamical, non-linear rescaling of time 
without changing the trajectories or equilibria from the standard replicator dynamics [1]. 

If we now focus on a homogenous population with type X individuals only, i.e. y = 0, in 
Eqs. (1), we have 

X = x{{a-\c)x - Ad) ■ (3) 

This equation has two equilibria, x = 0 and x = The first equilibrium x = 0, corresponding 
to extinction, is stable. The second equilibrium x = exists when a > Ac. However it is 
always unstable as x > 0 for x > and x < 0 for x < In this alternative model where 
interactions affect birth, a homogenous population either goes extinct or explodes, depending 
on the initial population size. Although, the deterministic equations Eqs. (1) and (2) appear to 
be reasonable, their ecological meaning remains unclear [2]. 


2 



2 Master equation and the diffusion approximation 


The microscopic process is a two dimensional Markov process in continuous time which can 
described by its master equation [3], 


dP{x,y,t) 

dt 


= T+_^P{x-l,y,t) + T+_^P{x,y-l,t) 
+ T~^iP{x + l,y,t) + T~^^P{x,y + l,t) 


(Tx P'^y T Tj; T ^P(^x^y^t^, 


( 4 ) 


where P{x,y,t) is the probability that there are x individuals of type X and y individuals of 
type Y at time t, and T is the transition rate of the population from one state to its neighbouring 
state. The subscript of T refers to the type whose density changes, and the superscript denotes 
whether its density increases by one or decreases by one. For example, T+ is rate that the 
number of type X increases from x to x + 1 and the number of type Y remains constant and 
r+ is rate that the number of type Y increases from y lo y + I and the number of type X 
remains constant. The transition rates can be deduced from the reaction rates and the number 
of individuals of those types involved in the corresponding reactions. 


T+ 

— ^x^xx X 


T+ 

y 

= \^yy y 


T~ 

X^ X y 

= \x^^x+ — + — 


T~ 

y 

, y^ X y 

( 5 ) 


Here, M is a scaling term which determines the frequency of competition compared to intrinsic 
growth. As we have shown in the main text, it controls the size of the system. When M is 
larger, the density of individuals in the deterministic equilibria is larger. 

To perform a diffusion approximation of the master equation, we scale the numbers x and 
y by M, X = x/M and y = y/M. The new variables x and y are approximately continuous for 
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sufficiently large M. We also rescale time ast = t/M. This leads to 


dP{x,y,t) 

dt 


- )^(^> t), 


( 6 ) 


Note when we rescale the time, we also need to rescale the transition rates accordingly. Thus 
in the same time unit, the transition rate from, for example, from state {x,y) to state {x + l,y), 
equals to the transition rate from the scaled state {x,y) and the scaled state (x + ■^,y)- This 
leads to 


r.+ 

X 

— MXx^xx X 


r+ 

y 

= M(Xy^yyy) 


r~ 

x 

, ~ 5?^ XU, 

= M(A.^ox + — + —) 


T- 

y 

- M{Xy^oy + ^ + ^)- 

( 7 ) 


Now we expand the transition rates and probability densities in Eq. (6) in a Taylor series at 
X or y, and we obtain 


Pix±^,y,t) 

Pix,y±^,t) 



TZ 1 

*+M 


y M 


TZ 1 
y~^ M 


~ dP{x,y,t) 1 d^Pix,y,t) 1 

’ ’ dx M dx'^ 2M‘^ ’ 


P{x,y,t) ± 


dP{x,y,t) 1 d^P{x,y,t) 1 


, dT+ 1 

'T'+ _ jZ 

* dx M 


dTz 

T~ + 

'Y- ' 


1 


dx M 

art 1 


TP - 
y dy M 


1 


dTz 

T~ H-^- 

y dy M 


dy M 
d^TP 1 
^ 5T2 2M2 


d^ 


2M2 ’ 


d'^Tz 


1 


5T2 2M2 

1 


52 


5 ^ 2M2 


d^Tz 


1 


dp 2M2 


( 8 ) 
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We denote P = P{x,y,t) and insert Eqs. (8) into Eq. (6), then we obtain 


dP 


+ 


+ 


+ 


T+ - 


erJ- 1 

dx M 


^ 9 T+ 1 

_ y__ 

^ y dy M 

^ dP- 1 

T~ H-- 

* dx M 



dy M 

\ 

T~ +T~ +T~ +T~ 




dP 1 
dxM~'^ 


dy M 


„ dP 1 

P H--h 

dx M 



dP 1 
dy M 


d^P 1 \ 

^2M2 ) 

d‘^P 1 \ 
~ff^2AP ) 

d'^P 1 \ 
WS2M^ ) 

d‘^P 1 \ 
1)^2M^ ) 


( 9 ) 


If we consider only the terms of the first order M~^ and the second order in Eq. (9), we 
obtain the Fokker-Planck equation 


dP 


+ 



d'^y \\ y y / J 



( 10 ) 


From Eqs. (7), we have -T~ = - \x^o)x --T~ = M{{Xy^yy - 

^y^o)y — ^ , T~ + T~ = M{{Xx^xx +^x^o)x + ^ + +T~ = M{{Xy^yy + 

Xy^Q)y + ^ + 1^). Putting these into Eq. (10), we can rewrite the Fokker-Planck equation as 


dP{x,y,t) 

dt 


dx * 




\^o)y cAf) 

19^/,, , xy\ 


( 11 ) 


The equivalent stochastic differential equations [4], which can often be handled in an numeri- 
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cally more efficient way, are 



where ^ is Gaussian white noise with mean 0 and variance 1. Note that the noise term vanishes 
when the population is close to extinction, but increases approximately linearly in the population 
density. 


3 The stability of the equilibria 


The deterministic rate equations for two types are 


X 

y 


X 




[ la: I y\ 

y\:--cM-dM) 


{13a) 

{13b) 


where x and y denote the numbers of individuals of type X and Y, and r refers to the same 
intrinsic growth rate for both types. For x = 0 and y = 0, we have four equilibria, Ei = (0,0), 
E 2 = {0,dMr), E 3 = {aMr,0), £’4 = ( ""tfSad '"fc-ad ^^e following, we perform a 

linear stability analysis of the four equilibria. The Jacobian matrix at the equilibrium {x*,y*) is 

r - ^ - Y- -YL 

' aM bM bM 

y* r — 

cM ' cM dM 

r 0 

(i) For £1 = (0,0), J{0,0) = . Thus, for any r > 0, the two eigenvalues are positive 

0 r 

and this equilibrium is unstable. 

r — tT 0 

with eigenvalues -r and -{d - b)r/b. If 

_ d rp _ iy> 

C ' 

d>b, the equilibrium is stable. Otherwise, it is unstable. 

—r —) 

(iii) For £2 = {oMr, 0), J{aMr,0) = 

0 r — 


r 

with eigenvalues -r and -{a - c)r/c. If 


(ii) For £2 = {OjdMr), J{0,dMr) = 



J{x\y*) = 
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a> c, the equilibrium is stable. Otherwise, it is unstable. 


(iv) For E 4 = Mr, the Jacobian matrix in this equilibrium is 




2 {b—d)c d{c—a) 

be—ad ^ be—ad ^ 

bd[c—a) 
c(bc—ad) ^ 

c{b—d) ac(b—d) 

be—ad b{be—ad) 

bd{e—a) b{e—a) 

e{be—ad) be—ad 


ae(b—d) 

~ b{bc-ad) ^ 
a{b—d) 2 b{e—a) 

be—ad ^ he—ad ^ 


( 15 ) 


Thus the eigenvalues A can be obtained given 

/ c{b - d) _ \ / b{c - a) _ \ _ ad{b - d){c - a) ^ 

v be — ad / \ be — ad ) (be — ad^ 

' —V—" '—V—" '——- 

A B c 


( 16 ) 


We can write Eq.(16) as - {A + B)X + AB -(7 = 0, thus A = ± \J + C. 

This equilibrium exists only if the numbers of individuals of both types are positive, i.e. 

> 0 and > 0- This yields two cases: b > d k e > a and b < d k, e < a, 

which both result in (7 > 0 and ^ + i? < 0. 

Thus, the first eigenvalue Ai = is always negative. The second 

eigenvalue A 2 = +~0 is negative if > {A-sf which can be 

simplified to AB - (7 > 0. From Eq. (16), we have 


AB-C 


e{b — d) b{e — a) ad{b — d){e — a) 
be — ad be — ad {be — ady 

{be — ad){b — d){e — a) 

{be — ady 
{b — d){e — a) 

{be — ad) 


( 17 ) 


Thus, W b > d k e > a, then AB - (7 > 0 and the equilibrium is stable; W b < d k e < a, 
the equilibrium is unstable. 


In general, although the equilibria in the deterministic limit are different from those in the 
replicator dynamics, their stability remains the same. 
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4 Ranking of equilibria in a coexistence game 


In a coexistence game where a < c and b > d, two types stably coexist with each other in 
the deterministic system. Here, we list all possible ranking rankings of the carrying capacities 
in two homogenous populations and the equilibrium densities in a heterogeneous population, 
see Tab. 1. One particularly interesting ranking is > Ry > > Kcox- a homogenous 

population of type X reaches higher numbers than type Y but in the mixed equilibrium type X 
is maintained at lower numbers than type Y. In stochastic processes this becomes particularly 
important because smaller carrying capacities result in a higher risk of extinction. For this 
particular ranking, a mutant Y has an increased chance to take over the entire population 
because in the mixed equilibrium X is outnumbered and hence is more likely to go extinct due 
to stochastic fluctuations. If this happens, the number of type Y individuals will increase until it 
reaches Ry. Over the course of this invasion, the total number of individuals changes from R^ 
to Rcox and ends at Ry. 

In co-existence games where b < a holds, the ranking is Ry < R"^ < Rqox and the total 
population size first increases and then decreases. In co-existence games where b > a holds, 
the ranking is Ry < Rcox < and hence the total population size continuously declines to the 
carrying capacity of type Y in isolation. This implies that evolution is not a simple process of ac¬ 
cumulating beneficial mutations, which have higher carrying capacities in isolation, but instead 
the adaptive process can favour invasion and fixation of mutations that are disadvantageous for 
the entire population [5]. Similar evolutionary patterns are apparent in the prisoner’s dilemma. 


Conditions Ranking 


a 

> 

d 

and 

ad{c 

-d) 

< 

bc{a - 

-d) 


RX 

> 

^c"ox > KV 

> 

Ry 

-‘'■COX 

a 

> 

d 

and 

ad{c 

-d) 

> 

bc{a - 

-d) 

and 

RX 

> 

Ry 

> R^ 

> 

-‘^cox 




b > c 

or 

ad{c 

- 

b) < bc{a 

-d) 







a 

> 

d 

and 

ad{c 

-b) 

> 

bc{a - 

-d) 


RX 

> 

Ry 

> R^ 

^ -‘'■COX 

> 

-‘'■cox 

a 

< 

d 

and 

ad{b 

- a) 

< 

bc{d - 

- a) 


Ry 

> 

K^cox > K- 

> 

RX 

-‘'■cox 

a 

< 

d 

and 

ad{b 

- a) 

> 

bc{d - 

-a) 

and 

Ry 

> 

RX 

V 

X 

> 

RX 

-‘'■COX 




b < c 

or 

ad{c 

- 

b) > bc{a 

-d) 







a 

< 

d 

and 

ad{c 

-b) 

< 

bc{a - 

-d)) 


Ry 

> 

RX 

V 

X 

> 

-‘'■cox 


Table 1: In a coexistence game (a < c and d < b), different conditions on the payoffs lead to 
certain rankings between the carrying capacities at homogenous and heterogenous equilibria. 
Here we list all possible rankings assuming identical intrinsic growth rates = ry. 
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5 Cyclic dynamics 


In our stochastic model, cyclic oscillation can be observed in a population with three types. 
Here, we show an example where the cycles spiral into an internal equilibrium according to 
the deterministic equations, but periodic oscillations persist under the corresponding stochastic 
process, see Fig. 1. 

According to the payoff matrix (see inset in Fig.1) in the deterministic dynamics, X can 
invade a homogenous Y population, Y can invade a homogenous Z population, and Z can 
invade a homogenous X population. The three types cycle into an internal equilibrium (see 
Fig. 1). However, demographic stochasticity drives the population away from the deterministic 
equilibrium and thus maintain the fluctuations over time. Note for small M, the population size 
is so small that demographic stochasticity will lead to the extinction of the whole population. 



Figure 1: Cyclic dynamics of three types under competitive Lotka-Volterra dynamics. The in¬ 
teractions are given by the payoff matrix in the inset. Saturated lines represent the deterministic 
dynamics, pale lines show one stochastic realisation (parameters Xx^xx = Xy^yy = X^^^z = 
0.6, = Xy^o = Xz-^o = 0.1, M = 10000, xq = 90, yQ = zq = 10). 
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6 Fixed population size 


Stochastic evolutionary dynamics under frequency dependent selection in populations of con¬ 
stant size, N, has recently attracted considerable interest, see e.g. [6, 7], The most popular ex¬ 
amples are based on birth-death models in discrete time. This results in a reduced set of micro¬ 
scopic interactions for two types, X+Y -> X +X, X+Y -^Y+Y. The deterministic rate equa¬ 
tions are x = (Xxy^xx-Kx^xy) X y and y = -x. In this case, x and y can be rescaled to indicate 
the frequencies due to constant population size. Constant reaction rates just lead to logistic 
growth of both types. Instead, to accommodate frequency dependent interactions, frequency 
dependent reaction rates are needed, e.g. Xxy^xx = fx/N and Xxy^yy = fy/N, where fx and 
fy represent the frequency dependent fitness. This yields x = x(l-x)(/x-/j/), which is just the 
replicator dynamics [8]. Alternatively, setting Xxy^xx = fx/ifX) and Xxy^yy = fy/{fN) where 
I = xfx + yfy denotes the average fitness of the population, leads to x = x(l - x){fx - fy)/] 
and recovers the adjusted replicator equation [9], which describes the deterministic limit of the 
frequency-dependent Moran process in discrete time [7]. Interpreting frequency dependent in¬ 
teractions under fixed population size based on individual reactions is somewhat problematic. 
First, it is unintuitive to define frequency dependent reaction rates on the microscopic level 
because each reaction should occur independently. Second, for an evolving population we ex¬ 
pect changes in the population size resulting by the changes of population composition rather 
than merely due to stochastic effects, which questions the basic tenet of models with constant 
population size. 
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